Flip dynamics in octagonal rhombus tiling sets 
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We investigate the properties of classical single flip dynamics in sets of two-dimensional random 
rhombus tilings. Single flips are local moves involving 3 tiles which sample the tiling sets via Monte 
Carlo Markov chains. We determine the ergodic times of these dynamical systems (at infinite 
temperature): they grow with the system size Nt like Cst. N^lnNr; these dynamics are rapidly 
mixing. We use an inherent symmetry of tiling sets and a powerful tool from probability theory, the 
coupling technique. We also point out the interesting occurrence of Gumbel distributions. 
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After the discovery of quasicrystals G[, quasiperiodic 
tilings as well as their randomized counterpart, ran- 
dom rhombus tilings ||, rapidly appeared to be suitable 
paradigmatic models for quasicrystalline alloys Q|. Si- 
multaneously, these systems also became an active topic 
in discrete mathematics (see || or jfj for examples). 
Fig. |l| displays a random tiling, which belongs to the 
class of plane tilings with octagonal symmetry. Beyond 
this case, plane tilings with larger symmetries including 
Penrose tilings Q and space tilings with icosahedral sym- 
metries were proposed to model every kind of quasicrys- 
tal. The present letter is devoted to dynamical proper- 
ties of random rhombus tilings in terms of local dynami- 
cal rules, the so-called phason-flips, which consist of local 
rearrangement of tiles (Fig. |l|). These dynamics are of in- 
terest for several reasons. On the one hand, it is more and 
more clear that phason-flips exist in real quasicrystals |7]] 
and can be modeled in a first approximation by tile-flips; 
they are a new source of atomic mobility, as compared 
to usual crystalline materials. In particular, they could 
carry their own contribution to self-diffusion even if 
the efficiency of such processes remains controversial || . 
They are also involved in some specific mechanical prop- 
erties of quasicrystals, such as plasticity via dislocation 
mobility jit)) . Therefore a complete understanding of flip 
dynamics is essential in quasicrystal physics. The present 
work is a first step in this direction. On the other hand, a 
lot of numerical work has been carried out to characterize 
statistical properties of tiling sets, a part of which was 
based on Monte Carlo techniques which rely on a faithful 
sampling of tiling sets (see |Tl| , ^2| ), So far, no systematic 
study of the relaxation times between two independent 
numerical measures has been accomplished, whereas it 
is an essential ingredient for a suitable control of error 
bars. However, there exist exact results in the simplest 
case of random rhombus tilings with hexagonal symme- 
try Jl3|-p^| and several estimates of relaxation times in 
larger symmetries, either numerical or in the approxi- 
mate frame of Langevin dynamics . 

Random rhombus tilings are made of rhombi of uni- 
tary side length. They are classified according to their 



global symmetries M. The simplest class of hexagonal 
tilings - made of 60° rhombi with 3 possible orienta- 
tions - has been widely explored {HS-flg. Tilings with 
octagonal symmetry are made of 6 different tiles (char- 
acterized by their shape and orientation): two squares 
and four 45° rhombi (Fig. |l|). Beyond these two cases, 
one can define tilings with higher symmetries (e.g. Pen- 
rose tilings (|] ) or of higher dimensions Q . For sake of 
technical simplicity, we focus on tilings filling a centrally 
symmetric polygon with integral side lengths (Fig. [l]). 
We are interested in the large size limit where the poly- 
gon becomes large keeping a fixed shape. Such tilings can 
schematically be seen as frozen near their boundary, and 
strain- free in their central region (Tt]] . Dynamics of fixed 
boundary tilings should be the manifestation of dynam- 
ics of their strain-free center, thus relating both bound- 
ary conditions. Here we suppose that all tilings have the 
same probability; we work at infinite temperature. 




FIG. 1. Examples of octagonal fixed boundary tiling and 
of elementary flip. We have also displayed (in gray) the de 
Bruijn lines of a family, among 4 families. They are lines of 
adjacent tiles sharing an edge with a fixed orientation. 

The set of all the tilings of such a region together with 
the flip dynamical rule define a discrete time Markov 
chain: at each step, a vertex of the tiling is uniformly 
chosen at random and if this vertex is surrounded by 3 
tiles in flippable configuration, then we flip it. Since sets 
of plane rhombus tilings are connected via elementary 
flips |TH |, this process can reach any tiling. It converges 
toward the uniform equilibrium distribution, since it sat- 
isfies detailed balance. All the difficulty is to character- 
ize how many flips one needs to get close to equilibrium. 
Generally speaking, let us consider a Markov chain on 
a finite configuration space L, which converges toward a 



1 



stationary distribution 7T. Let xq be any initial configura- 
tion and P(x, t\xo, 0) be the probability that the process 
has reached the configuration x after t steps. Then 

A(t,x ) = l/2 Y,\ p (x,t\*o,0)-*(x)\ (1) 

usually measures the distance between both distribu- 
tions [fl9|| . Given e > we define the ergodic or mixing 
time r(e) so that whatever xq, after r(e) steps, one is 
sure to stay within distance e of equilibrium: 

r(e) = maxmhi{£o/Vi > t ,A(t,x ) < e} . (2) 

x t 

In this letter, we prove that for a tiling of Nt tiles 

t(s) < Cst. (N T )»\nN T \n(l/e). (3) 

More precisely, we establish, with the help of reduced nu- 
merical work, that such a bound holds for v = 3, then we 
argue that v — 2 should be the correct exponent. What- 
ever this exponent, this proves that flip dynamics are 
rapidly mixing at infinite temperature. As for the ln(A^) 
correction in (||) , as discussed in jl5| , it is a feature of our 
choice of distance A(i, xq): an Euclidean norm would not 
display this correction but it is less natural in the context 
of measure of probability distributions convergence. 

The coupling technique |n| has been successfully ap- 
plied to estimate mixing times of several systems, such as 
hexagonal tilings (l^ [H| . It relies on the surprising idea 
that following the dynamics of couples of configurations 
instead of a single one might provide the properties of the 
original dynamics on single configurations. A coupling is 
a Markov chain on L x L; couples of configurations are 
updated simultaneously and are strongly correlated, but 
each configuration, viewed in isolation, performs transi- 
tions of the original Markov chain. Moreover, the coupled 
process is designed so that when both configurations hap- 
pen to be identical, then they follow the same evolution 
and remain identical forever. Then the central idea of the 
technique is that the average time the two configurations 
need to couple (or to coalesce) provides a good upper 
bound on the original mixing time r(e): given an initial 
couple (xo,yo) at time t = 0, define the coalescence time 
T(xo,yo) as the minimum time that both configurations 
need to coalesce, and the coupling time as 

T= max <T(x ,y )>, (4) 

{x ,y )eLxL 

where the last mean is taken over realizations of the cou- 
pled Markov chain. The following theorem, central in 
the coupling approach, provides an upper bound for the 
ergodic times of the original {not coupled) process fL9|| : 

r(e) < Teln(l/e) + 1 ~ Teln(l/e), (5) 

where e = exp(l). If the configuration set L can be en- 
dowed with a partial order relation y with unique min- 
imum and maximum elements and 1, the implementa- 
tion of the technique is highly facilitated, provided the 



coupled dynamics is monotonous, i.e. if x(t) y y(t), then 
x(t +1) h y{t + 1): let (xo, ?/o ) be any initial couple such 
that 1 y xq y yo y 0. After any number of steps, the 
4 configurations remain in this order. When the iterates 
of and 1 have coalesced, the iterates of xo and yo also 
have. Thus T(x ,y ) < T(6, 1) and T = < T(6, 1) >. 

A convenient representation of random rhombus tilings 
was introduced by de Bruijn p0[ . It consists of following 
in a tiling lines made of adjacent tiles sharing an edge 
with a fixed orientation (Fig. [l]) . The set of lines associ- 
ated with an orientation is called a de Bruijn family. In 
an octagonal tiling, there are 4 families. When removing 
a family from an octagonal tiling, one gets an hexagonal 
tiling. Conversely, this remark enables one to propose 
a convenient construction of tilings J2l],[l2|;|l8| : directed 
paths are chosen on an hexagonal tiling, called the base 
tiling. They are represented by dark lines in Fig. |[ They 
go from left to right without crossing (but they can have 
contacts) . When they are "opened" following a new edge 
orientation, they generate de Bruijn lines of the fourth 
family. In this one-to-one representation, a tiling flip in- 
volving tiles of the fourth family becomes a path flip: the 
path jumps from one side of a tile to the opposite side. 




FIG. 2. Example of 2-line coupling in the directed path 
representation. Dark lines represent a configuration and dot- 
ted lines the second one. Both configurations form a couple. 

Defining couplings on the whole octagonal tiling sets 
seems to be an infeasible task. We instead use the paths- 
on-tiling point of view to decompose the configuration 
space into smaller subsets where couplings can be defined: 
let J a denote the set of tilings which have the same base 
hexagonal tiling a, called "fibers" JTsfl . L is a disjoint 
union of fibers. The only possible flips inside J a are those 
which involve the fourth de Bruijn family. Note that we 
can construct four such fibrations, one for each family. 

Now, let Ai denote the symmetric transition matrix as- 
sociated with the Markov chain on the whole set L: given 
two configurations x and y, the matrix entry A4(x, y) is 
equal to the transition probability P(x, t + l\y, t). In the 
same way, we define the symmetric transition matrices 
M.i associated with the Markov chains where only flips 
involving the i-th de Bruijn family are allowed. Since 
fibers have been disconnected, M.i is block-diagonal. The 
following result interconnects the four fibrations: 

M = (M 1 +M 2 +M 3 + M 4 )/3 -Id/3, (6) 
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where Id is the identity. Indeed, each coefficient M(x, y) 
appears in all four matrices M.i but one, since the corre- 
sponding flip involves 3 de Bruijn lines. 

Now we implement the above coupling technique on 
each fiber. To begin with, we suppose that there is only 
one line in the flipping family, denoted by £. As in ref- 



erence 13 1, in order to have a monotonous coupling, we 
slightly modify the Markov chain: at each step, we choose 
uniformly at random an internal vertex of £, the n-th one 
starting form the left, and a number r G {0, 1}. If this 
vertex is fllppable upward (resp. downward) and r — 
(resp. r = 1), then we flip it. Note that this Markov 
chain has a time unit different from the original one. 

We now define the order relation (V): given two lines £\ 
and £2, l\ >z £2 if £\ is entirely above £2- The maximum 
(resp. minimum) configuration clearly lies on the top 
(resp. bottom) boundary of the hexagonal domain. If 
the two flips on £\ and I2 occur with same n and r and 
if l\ y £2 then their images satisfy the same relation. 
Indeed, as in reference thanks to the introduction 
of r, if a flip could bring l\ below £2, then the same flip 
would also apply to £2, thus preserving the order between 
lines: the coupling is monotonous. 

In the general case with p non-intersecting lines in each 
configuration (Fig. ||) , let us denote by £^ , j = 1, . . . , p, 
the p lines of each configuration 7,. Then 71 >; 72 if for 



y £2 ■ The configuration 7 is maximum (resp. 



anyj,^ 

minimum) when each of its lines is maximum (resp. min- 
imum). At each step, the index j of the line to be flipped 
is chosen between 1 and p, the same j for both 7^. 

To begin with, we numerically study the diagonal case, 
where the 4 sides of the octagonal tilings are equal to k. 
For a given base tiling 7^, we run a number m of cou- 
plings until they coalesce, and then estimate the coupling 
time T(T a ). We then make a second average on M differ- 
ent tilings, in order to get the time T averaged on tilings 
T a . We also keep track of the standard deviation AT. 
From our numerical data (see Fig. ||), we draw the fol- 
lowing conclusions: AT/T decreases toward a constant 
(~ 0.07) as k — > 00, which means that the average cou- 
pling time T(T a ) goes on depending on the base tiling 
T a at the large size limit. However, most T(T a ) are of 
order T, and the mixing times t(e) on most fibers are 
controlled by T. Nevertheless, the effect of few "slower" 
fibers will deserve a detailed discussion below. Moreover, 
the measures of T are compatible with a k 4 In k behavior 
(Fig. |^, inset). In particular, this fit with logarithmic 
corrections is much better than a simple power-law fit. 
This result is consistent with known results in the case of 
hexagonal tilings |K^ , |l5| ], where T also grows like A; 4 lnfc. 

We also have explored coupling times on fibers in non- 
diagonal cases and our conclusions remain identical. As 
a consequence, couplings in fibers behave like couplings 
in hexagonal tiling problems, up to different numerical 
prefactors: the dynamics on each fiber is rapidly mixing. 
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FIG. 3. A distribution of coupling times T in the case of a 
diagonal base tiling (k = 10) and p = 1, fitted by a Gumbel 
distribution (continuous curve; see eq. ([]))• Inset: in the 
diagonal case (k — p), numerical estimates of T in function of 
k 4 \nk (circles), up to k = 15, and linear fit. Error bars are 
smaller than the size of symbols. The slope is 25.51 ± 0.05. 




FIG. 4. The lattice of tilings filling an octagon of sides 1,1,1 
and 2. The edges represent possible flips. Two fibrations 
among four are represented (continuous and dotted lines). 

Now we return to the dynamics on the whole set of 
tilings. Examining Fig. ^, one remarks that two fibra- 
tions are to a certain extent "transverse" : in octagonal 
tiling sets, one can connect any two tilings using flips of 
only two fibers p2| ]; if the dynamics is rapidly mixing 
in each fiber, the combination of dynamics on two (and 
even four) fibrations will certainly also be rapidly mixing. 
We now establish properly this point: it is common in the 
field of Markov processes to relate rates of convergence to 
spectra of transition matrices. Generally speaking, given 
a transition matrix M, 1 is always the largest eigenvalue 
in modulus, and the difference g(M) between 1 and the 
second largest eigenvalue is called the first gap of M. 
Then r(e) ~ ln(l/ 'e)/g(M) for small e @. By simple 
arguments from linear algebra and Euclidean geometry 
(since matrices are symmetric), the following central gap 
relation can be established J22|, based upon @ and the 
smallness of the intersection of two fibers: 

g(M)>M(g(Mi)). (7) 
This result restores the symmetry lost in the fibration 
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process. It implies that the mixing time r(e) on L is 
smaller than the mixing time on the slowest fibration de- 
rived from (||). Hence r(e) is smaller than the mixing 
time on the slowest fiber. We have seen that the aver- 
age coupling time T grows like fc 4 lnfc in diagonal cases. 
Since the number of tiles is Nt — 6k 2 , T ~ k±N t Ih(Nt), 
where K4 = 1.189 ± 0.003 in the original time unit. How- 
ever, the coupling time T depends on the base tiling T a 
and the distribution of times T in a given fibration have a 
certain width around T (Fig. ||). And even if the typical 
values of coupling times are of the order of magnitude 
of the previous average value, this does not exclude the 
existence of rare slow fibers in the upper tails of these 
distributions. However, we recall that coupling times T 
are maxima of coalescence time distributions (^) . There- 
fore the expected shape of the distribution of T(T a ) ought 
to be sought in the specific class of extreme- value distri- 
butions, namely Gumbel distributions p3|: consider N 
independent identical random variables T a , whose prob- 
ability densities decay rapidly at large T: 

p(T)^exp(-C 2 T^), (8) 

where C\,C2,/3 > 0. If T max = max a T a , then at large 
N, the probability density of T max satisfies 

p(u) = exp( — u — exp(— u)), (9) 

where u = (T max — To) /ST is a suitably rescaled variable. 

Now, even if coalescence times are not strictly speak- 
ing independent variables p^| , our numerical distribu- 
tions appear to be well fitted by this kind of distribu- 
tion (Fig. This result provides the large T behavior 
of coupling time distributions: p(T max ) ~ exp(— u) ~ 
exp(— T max /^T), and therefore an estimation of the 
largest coupling time. Let us focus on diagonal cases: 
To as well as ST behave like k 4 In k. But for a fibration 
i, there are Ni base tilings T a . Therefore if T* is the 
largest coupling time on all tilings T a , it is estimated by 
Ni exp(— T*/5T) w 1. Now Ni grows exponentially with 
the number of tiles @: In A^ w Csti k 2 . Thus 

T* « Cst 2 k 6 In k w Cst 3 N T In Nt (10) 

and t(e) < Cst^ N T In Ay ln(l/e). However these ex- 
treme values should not be relevant: because of the ex- 
ponential decay of p{u) , slow fibers are rare and can be 
bypassed via rapid ones. More precisely, perturbation 
theory arguments suggest that rare slow fibers can be 
seen as a small perturbation of an otherwise rapid transi- 
tion matrix and have a vanishing influence on its spectral 
gap p2]: they do not significantly slow rapid dynamics 
and the typical coupling time T drives the dynamics on 
fibers, leading to v = 2 in (|J) and Cst. w e K4. 

What does this analysis become in the case of larger 
symmetry tilings or of higher dimensional tilings, such 
as icosahedral ones? Plane rhombus tilings with 2Z?-fold 



symmetry could be addressed by our approach without 
significant technical complication, leading to laws similar 
to (0), up to different prefactors Cst. p2| . As for higher 
dimensional tilings, the fibration process remains valid, 
but the connectivity of fibers is not established jl8j , mak- 
ing impossible a naive generalization of this approach. 

I would like to thank M. Latapy, K. Frahm, R. Mosseri 
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